MODULE fem_main

   USE fem_global_types
   USE fem_global_consts
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fem_initial,                    &
             fem_final,                      &
             fem_get_boundary_phi,           &
             fem_solve_laplace_equation

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fem_initial(mode,ndim,lambda)

      USE fem_scheme_module, ONLY: fem_init_scheme,       &
                                   fem_get_scheme
      USE fem_A_matrix,      ONLY: fem_init_A_matrix

      IMPLICIT NONE
      CHARACTER(*),  INTENT(IN) :: mode
      INTEGER(INTK), INTENT(IN) :: ndim(3)
      REAL(REALK),   INTENT(IN) :: lambda(3)

      TYPE(fem_scheme), POINTER :: scheme

      WRITE(IO1,*)
      WRITE(IO1,*) "-----------------------------------------"
      WRITE(IO1,*) "|         FEM initialization            |"
      WRITE(IO1,*) "-----------------------------------------"
      WRITE(IO1,*)
      CALL fem_init_scheme(0)
      CALL fem_get_scheme(scheme)
      CALL fem_setup_geometry(scheme,ndim,lambda)
      CALL fem_init_A_matrix(scheme)

   END SUBROUTINE fem_initial

!-------------------------------------------------------------------------------

   SUBROUTINE fem_final

      USE fem_A_matrix, ONLY: fem_free_A_matrix
      IMPLICIT NONE
      CALL fem_free_A_matrix

   END SUBROUTINE fem_final

!-----------------------------------------------------------------------------

   SUBROUTINE fem_setup_geometry(scheme,ndim,lambda)

      IMPLICIT NONE
      TYPE(fem_scheme), INTENT(INOUT) :: scheme
      INTEGER(INTK),    INTENT(IN)    :: ndim(3)
      REAL(REALK),      INTENT(IN)    :: lambda(3)

      ! Single element properties
      scheme%FEFNNod1D = scheme%FEFNthShp + 1
      scheme%FEFNNod3D = scheme%FEFNNod1D ** 3
      scheme%FEFNodEdg(:) = lambda(:)

      ! Whole solution space
      scheme%FEFNPt1D(:) = ndim(:)
      scheme%FEFNPt3D = PRODUCT(ndim(1:3))
      scheme%FEFNElm1D(:) = (ndim(:)-1) / scheme%FEFNthShp
      scheme%FEFNElm3D = PRODUCT(scheme%FEFNElm1D(1:3))
      ! Interior points only
      scheme%FEFNNB3D = PRODUCT(ndim(1:3) - 2)

      WRITE(IO1,*) "scheme%FEFNthShp    =",   scheme%FEFNthShp
      WRITE(IO1,*) "scheme%FEFNNod1D    =",   scheme%FEFNNod1D
      WRITE(IO1,*) "scheme%FEFNNod3D    =",   scheme%FEFNNod3D 
      WRITE(IO1,*) "scheme%FEFNodEdg(:) =",   scheme%FEFNodEdg(:)
      WRITE(IO1,*) "scheme%FEFNElm1D(:) =",scheme%FEFNElm1D(:)
      WRITE(IO1,*) "scheme%FEFNElm3D    =",   scheme%FEFNElm3D
      WRITE(IO1,*) "scheme%FEFNPt1D(:)  =", scheme%FEFNPt1D(:)
      WRITE(IO1,*) "scheme%FEFNPt3D     =",    scheme%FEFNPt3D
      WRITE(IO1,*) "scheme%FEFNNB3D     =",    scheme%FEFNNB3D

      ! Error checks
      IF ( SUM( MOD(ndim(:)-1, scheme%FEFNthShp) ) /= 0) THEN
         STOP 'boundary points not exact multiple of FE nodes'
      END IF
      
   END SUBROUTINE fem_setup_geometry

!-------------------------------------------------------------------------------

   SUBROUTINE fem_get_boundary_phi(bc_file,phi)

      IMPLICIT NONE
      CHARACTER(*), INTENT(IN)  :: bc_file
      REAL(REALK),  INTENT(OUT) :: phi(:,:,:)

      INTEGER(INTK) :: i, nx,ny,nz, npoints
      REAL(REALK) :: tmp
      LOGICAL :: fail

      OPEN (UNIT=IO2,FILE=bc_file,STATUS='OLD',     &
            ACTION='READ',FORM='UNFORMATTED')
      REWIND (IO2)

      READ(IO2) npoints
      fail = .FALSE.
      DO i = 1, npoints
         READ(IO2) nx, ny, nz, tmp
         IF (nx > SIZE(phi,1)) fail = .TRUE.
         IF (ny > SIZE(phi,2)) fail = .TRUE.
         IF (nz > SIZE(phi,3)) fail = .TRUE.
         phi(nx,ny,nz) = tmp
      END DO

      IF (fail) STOP 'inconsistency in boundary points!'
      CLOSE(UNIT=IO2,STATUS='KEEP')

   END SUBROUTINE fem_get_boundary_phi

!-------------------------------------------------------------------------------

   SUBROUTINE fem_solve_laplace_equation(bc_file,phi)

      USE fem_scheme_module, ONLY: fem_get_scheme
      USE fem_linear_solver, ONLY: fem_init_linear_equations,    &
                                   fem_free_linear_equations,    &
                                   fem_solve_linear_equations

      IMPLICIT NONE
      CHARACTER(*), INTENT(IN)  :: bc_file
      REAL(REALK),  INTENT(OUT) :: phi(:,:,:)

      TYPE(fem_scheme), POINTER :: scheme

      CALL fem_get_scheme(scheme)

      ! Initialize phi from boundary condition
      phi(:,:,:) = zero
      WRITE(IO1,*) "Call fem_get_boundary_phi..."
      CALL fem_get_boundary_phi(bc_file,phi)

      ! Initialize and solve Laplace equation for phi
      WRITE(IO1,*) "Call fem_init_linear_equations..."
      CALL fem_init_linear_equations(scheme,phi)
      WRITE(IO1,*) "Call fem_solve_linear_equations..."
      CALL fem_solve_linear_equations(scheme,phi)
      CALL fem_free_linear_equations

   END SUBROUTINE fem_solve_laplace_equation

!-------------------------------------------------------------------------------

END MODULE fem_main
